% 求解函数最小值 $F(x) = \sum\limits^D_{i=1} x_i^2$
clear
close all
clc
 
SearchAgents_no = 30 ; %　种群规模
dim = 10 ; % 粒子维度
Max_iter = 1000 ; %　迭代次数

% 约束上下界（本题中是x的取值范围）
ub = 5 ;  % 约束上界
lb = -5 ; % 约束下界

% 学习因子
c1 = 1.5 ; %　学习因子1
c2 = 1.5 ; %　学习因子2

% 惯性权重
w = 0.8 ;
% 粒子最大速度
vmax = 3 ; 

% 取随机数进行初始化（取种群规模数量的随机数，代入公式作为位置和速度）
pos = lb + rand(SearchAgents_no,dim).*(ub-lb) ; % 初始化粒子群的位置
v = - vmax +2*vmax* rand(SearchAgents_no,dim) ; % 初始化粒子群的速度

% 初始化每个历史最优粒子

% 认为pos是当前最优，并存储到pBest中
pBest = pos ; 
% 用一维向量存储整个种群当前状态的最优值（适应值）
pbestfit = zeros(SearchAgents_no,1);
for i = 1:SearchAgents_no
    pbestfit(i) = sum(pos(i,:).^2) ; 
end

% 初始化全局历史最优粒子

% 找到种群向量中的最小值（最佳粒子），存储在gBest中
[gBestfit,index] = min(pbestfit) ;
gBest = pos(index,:) ;
Convergence_curve = zeros(Max_iter,1);
 
% 开始粒子群算法

for t=1:Max_iter
    for i=1:SearchAgents_no
        % 更新个体的位置和速度
        v(i,:) = w*v(i,:)+c1*rand*(pBest(i,:)-pos(i,:))+c2*rand*(gBest-pos(i,:)) ;
        pos(i,:) = pos(i,:)+v(i,:) ;
        % 边界处理
        v(i,:) = min(v(i,:), vmax);
        v(i,:) = max(v(i,:), -vmax);
        pos(i,:) = min(pos(i,:), ub);
        pos(i,:) = max(pos(i,:), lb);

        % 更新个体最优
        
        % 计算当前位置（pos）的适应值（f1），并与原适应值（pbestfit）进行比较，判断是否更优
        f1 = sum(pos(i,:).^2);
        if f1 < pbestfit(i)    
           pBest(i,:) = pos(i,:) ;
           pbestfit(i) = f1;
        end
        % 更新全局最优

        % 代码逻辑同上
       if pbestfit(i) < gBestfit
            gBest = pBest(i,:) ;
            gBestfit = pbestfit(i) ;
       end
    end
    % 每代最优解对应的目标函数值
    Convergence_curve(t) = gBestfit; 
    disp(['Iteration = ' num2str(t)  ', Evaluations = ' num2str(gBestfit)]);
end
 
% 画图

figure('unit','normalize','Position',[0.3,0.35,0.4,0.35],'color',[1 1 1],'toolbar','none')
subplot(1,2,1);
x = -5:0.1:5;y=x;
L=length(x);
f=zeros(L,L);
for i=1:L
    for j=1:L
       f(i,j) = x(i)^2+y(j)^2;
    end
end
surfc(x,y,f,'LineStyle','none');
xlabel('x_1');
ylabel('x_2');
zlabel('F')
title('Objective space')
 
subplot(1,2,2);
semilogy(Convergence_curve,'Color','r','linewidth',1.5)
title('Convergence_curve')
xlabel('Iteration');
ylabel('Best score obtained so far');
 
axis tight
grid on
box on
legend('PSO')
display(['The best solution obtained by PSO is : ', num2str(gBest)]);
display(['The best optimal value of the objective funciton found by PSO is : ', num2str(gBestfit)]);
 
        
 
 
 
